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ABSTRACT 


This thesis presents a computer program which will 
solve for the lowest natural vibrational frequency of an 
arbitrarily shaped three dimensional object. The program 
was developed using a finite element numerical model which 
may employ either linear or quadratic isoparametric elements. 


Necessary instructions and guides for the user are included. 
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Pee hn RODUCT LON 


This thesis presents and discusses a computer program, 
called NAT1-3D, which employs the finite element method to 
set up and solve the elastodynamic equations for the free 
Vibration of a three-dimensional, linearly elastic body 86 
as to emerge with an accurate approximation of the funda- 
mental frequency of free harmonic vibration. 

The investigation was originally motivated by the 
interesting, important, and unsolved problem of the effect 
of foundation compliance and foundation inertia upon the 
funcamental natural frequency of a short elastic cantilever 
meam. lt 1s clear that the usual assumption that the root 
of such a beam is rigidly fixed does not adequately represent 
miemract that mo foundation 1s truly rigid. Calculations 
based upon ideal root fixity thus necessarily result in 
estimates of the fundamental vibration frequency that are 
too high. There seems to be no reasonably reliable way of 
Baieimating this arraer Rei. lle. 

During the course of the work involved in adapting 
existing computer software to the problem indicated in the 
last paragraph (in particular adding a dynamic capability 
which such software did not previously possess) it became 
more and more apparent that it would be most useful to 


focus attention upon the development of a general purpose 


References will usually be cited in this manner; 
seecepage 6/7 Tor list. 
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program for the evaluation of the fundamental frequency cof 
three dimensional elastic bodies of any shape. Because of 
the usual and not unexpected problems of testing and debugging 
a very complicated digital computer program, time was not 
available for the study of the root fixity problem. Further- 
more, and more importantly, for this problem or any problem 
in which truly complicated geometry is involved, digital 
computer programs which rely solely on in-core storage, as 
does the program developed herein as well as most of its 
predecessors, are capable of supplying results of marginal 
accuracy so that parametric studies or studies involving 
systematic changes in geometry are apt to have limited 
mugeess. (his is because in-core storage limitations do not 
permit subdividing the body into sufficiently small parts 
to permit obtaining really accurate results. Accordingly, 
one of the recommendations made herein is that attention be 
focussed on modifying this program so as to permit efficient 
out-of-core storage and manipulation. The present study does 
demonstrate that the procedural sequence adoptea and imple- 
mented herein leads to reliable evaluations if the geometry 
of the vibrating object is sufficiently simple that it may 
be represented by a reasonably small number of finite 
elements. 

Although the immediately following sections (II, III, IV, 
Wemrand VI} discuss theory, sources, and details of the program 


gevelopec herein, the principal thrust and value of this 








thesis is felt to be the detailed instruction provided for 
the employment of the program NAT1-3D by a new user for the 
actual evaluation of the fundamental frequency of a three- 


dimensional elastic body in which he may be interested. 





ri, “Wwe ER UMetc ELEMENT METHOD 


ee 


In recent years a procedure, called the Finite Element 
Method, has been developed and widely used for a variety of 
engineering problems involving the analysis of bodies having 
complicated shapes. In the analysis of the stresses and 
deformations of an elastic body, using this method (which is 
generally abbreviated as FEM), the body is subdivided into 
a (large) number of geometrical subdivisions called finite 
elements. The adjoining surfaces of these elements are 
delineated by points called nodes. In a three-dimensional 
problem each node may have three displacements. Displace- 
momes. dnd, bY a process of differentiation, strains at 
points interior to such a finite element are expressable in 
terms of nodal displacements by use of so-called shape 
functions which, roughly speaking, are polynomials in the 
nodal displacements. The laws of elasticity are used to 
relate nodal displacements of a finite element to the forces 
applied at its nodes and the salavioasiip is conveniently 
expressed in the form of an element stiffness matrix. 

An assembly algorithm is used to combine such element 
stiffness matrices into a so-called global siffness matrix, 
merce original entire body. This global siffness matrix, 
hereinafter designated as [BGK], relates nodal displacements 


for the entire body to nodal forces which may be applied to 
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it. If there are a large number of finite elements in the 
body and if each of them has itself a substantial number of 
nodes, the total number of displacements, three for each node, 
may be very large. 

Im Clastostatic problems, the nodal forces, 1.e., the 
force components applied to the nodes, are arranged in the 
form of a vector and the unknown nodal displacement components 
are similarly arranged, so that the relationship can be 


expressed in the matrix equation 
[BGK] {v} = {Ff} 


Memenemsize of the vectors {v} and {f} were not great, any 
of a number of procedures could be used to solve this matrix 
equation for the unknown vector {v}. However, in practice, 
the size usually is quite large, and severe problems of 
storage in an electronic computer have resulted in the 
development of quite subtle procedures for storing the 
elements of [BGK] and of {f}. 

The most significant features of the matrix [BGK] are 
that it is positive definite, symmetric, and of banded 
Structure. The last means that there is an upper right and 
a lower left triangle which consists entirely of zero elements 
which need not actually be stored if their presence is other- 
wise suitably noted. These properties permit storing the 
mee) matrix amefar less Space than if it were to be stored 


tin the form usually used to exhibit square matrices. 
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Principle features of any FEM computer program are the 
Peeation and compact storage of the elements of [BGK], and 
a so-called "equation solver" which efficiently solves the 
linear algebraic problem indicated by the equation above. 
Software to accomplish this has been developed and tested 
andis readily available. One of the problems facing a person 
who undertakes to solve an engineering problem by FEM is to 
select among several competing software programs. JRISOP 
[Ref. 3] and FIELD [Ref. 4] are the names of programs 
developed by Professor G. Cantin, LCDR E.A. Leonidas and 
LCDR G.T. Lew. These programs have been successfully 
employed at NPS for a variety of investigations. It was 
decided to adapt and employ appropriate sections of these 
Baodrans for the purposes of this thesis. 

TRISOP and similar programs require a massive input to 
describe the geometry of the body and the manner in which it 
is subdivided, to introduce the physical properties of the 
materials, to represent the manner in which the body is 
@onstrained, and to account for the applied loadings. For 
example, in a typical problem treated in the course of this 
thesis investigation 250 data cards would be required as 
input. Since the labor of generating the data for and 
actually punching these cards is so great and since the 
likelihood of making at least one crippling error is almost 
certain, software programs called "mesh generators" have 


been developed which accept a minimum of input data and 
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MmcsrrueLions andeawnich prodtucesdata cards suitable for 
subsequent use in an FEM program. Such a mesh generator, 
called TRIMEG [Ref. 2], developed by LT J.R. Adamek, for use 
with TRISOP has been used and tested at NPS, and has provided 
a very useful, almost essential adjunct to the FEM program 
ised for this thesis. 

TRISOP does not have a dynamic capability. As will be 
indicated in the next section hereof, one of the important 
features of a dynamics FEM program for use in determining 
natural frequencies is an appropriate mass matrix. Without 
going into detail, it may be stated that the theory of the 
FEM provides for the construction of so-called consistent 
mass matrices for the individual finite elements and their 
assembly into a global mass mass matrix [BGM]. The computa- 
tional details are almost exactly the same as for the 
construction of element stiffness matrices and their assembly 
into [BGK]. The global mass matrix has essentially the same 
properties as the global stiffness matrix and requires the 
Same storage. 

Much thought was given to the problems of storing such 
a consistent mass matrix. However, studies by others 
[Ref. 7] have indicated that when evaluating the fundamental! 
harmonic frequency of free vibrations of an elastic object, 
use of a lumped-mass form of the mass matrix gives results 
as good as or better than those obtainable by using a con- 


Sistent mass matrix. The reasons for this are not clear, 


lee 





but the advantages are. A lumped-mass mass matrix is 
diagonal in form and can be stored as a vector; briefly, 
using the lumped-mass form results in a dynamics problem 
requiring very little more storage than an elastostatic 
problem requires. 

Accordingly it was decided to employ the lumped-mass 
form and the question arose of how optimally to obtain such 
a mass matrix. The procedure which finally evolved was to 
construct the element consistent mass matrix and to diago- 
nalize it before assembly into the global mass matrix. 

The word diagnaiize in this context is not intended to 
mean the mathematical process of diagonalization by means of 
solving the complete eigenvalue problem for the matrix M. 
Instead, it merely means replacing the consistent (element) 
mass matrix by a diagonal matrix corresponding to an 
appropriate division of the mass of the element into submasses 
concentrated at the nodal points. 

The assembly was dane in the core storage space reserved 
for assembling (and later operating upon) the global stiffness 
matrix, and the assembly algorithm is exactly the same for 
both matrices. The assembled mass matrix is diagonal and 
its nontrivial elements are quickly stored elsewhere in the 
fern Of a vector. 

The details of diagonalizing the element mass matrix 
deserve a brief description, explanation, and plausible 


justification. Aside from the fact that it is obviously 





desirable to employ a procedure which involves a simple 
algorithm and results in the same diagonalized element mass 
matrix regardless of renumbering of nodes, a reasonable and 
plausible distribution of the total element mass into mass 
points at the nodes is also required. This matter is 
discussed in somewhat greater detail on page 29. 
Accordingly appropriate segments of the programs TRISOP 
and FIELD were selected and modified so as to incorporate 
the construction of a lumped-mass mass matrix as outlined 
in the preceding paragraphs and so as also to incorporate 
the operations of LDLT decomposition, matrix (Stodola) 
iteration, convergence testing, and Rayleigh quotient 
evaluation discussed in the following section hereof. The 


resulting program is called NAT1-3D. 





TIT. METHOD OF SOLUTION 


— 





With the use of the [K] and [M] matrices provided by 
the finite element formation, as discussed previously, the 


free vibration problem can be formulated as 
[M] {v} + [K] {v} = {0} (1) 


The remarks which follow immediately represent standard 
procedure in vibration theory and are intended merely to 
recall to the reader the operational sequence which leads to 
the desired result. We are concerned with the lowest 
frequency harmonic vibrations which the system may execute. 
By assuming that the vector {v} of time-varying displacements 


has the form 
ies Cos ( wit) 


(where {v*} is a vector of constant quantities representing 
the amplitudes of the motions of the nodal points) and 
substituting this equation into equation (1) above, the 
cosine term may be cancelled out, and, upon dispensing with 
Meeeacterisk and henceforth regarding {v} as a vector of 


constants, we arrive at the form 


[Kk] {v} = w* [M] {v} 





It may be shown that if an arbitrary vector fs us 
introduced on the right side and this equation is solved 
for a new vector, 1.e., if we solve for re in the 


equation 
[Kk] {v.41} = w°[M] (v,} 


then ive 44} is a better approximation to the theoretically 
correct vector for the lowest frequency than was {v.}. 
Furthermore, if the vectors are suitably and similarly 
normalized, the scalar constant of normalization approaches 
more and more closely to the correct value of the square of 
the lowest frequency. This process of matrix iteration is 
frequently called Soon cea gee Paleaion after the Swiss engineer 
who introduced the method. 

Although the scalar constant of normalization does 
provide an increasingly good evaluation of the desired value 
Ot - at any stage of the iteration a more accurate estimate 


may be obtained from the Rayleigh quotient 
we = Q= tvd! Ck] tv / vd! CM] tv}. 


Performing this operation after a reasonable number of 
iterations provides the best currently available value for 
ws and, presuming that the value so obtained is "close" to 
the values previously obtained during the Stodola iteration, 
provides a check on the accuracy and the validity of the 


entire computational process. 


17 





With systems having only a few degrees of freedom it 
is customary to solve for Oy Pe by obtaining the inverse of 
the matrix [—K]. However, the desirable property of the 
[K] matrix as obtained by the finite element method, namely 
its banded character which permits compact storage, is not 
possessed by the inverse ree furthermore the operation 
of obtaining an inverse 1S computationally extravagant. 
Accordingly, it is very advantageous to represent [K] in the 
so-called LDLT form, i.e., to decompose [K] into the matrix 


product 
[kK] = [L] fo] 1]! 


where [L] is a lower triangular matrix having unit diagonal 
elements (which need not be stored) and [D] is a diagonal 
matrix. This decomposition is a standard operation which is 
noemaifficult to perform and the matrices [L] and [D] may 

be stored in the space occupied by [K]. 

The normalization which is used is the so-called 
Euclidean normalization for which the unnormalized vector 
{v} is divided throughout by the square root of the sum of 
the square of its elements. 

The iteration proceeds as follows. Having the normalized 
vector {v.} (either as a result of having made a previous 
iteration or, the first time through, as the result of a 


guess or enlightened estimate) the vector 
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is obtained simply by matrix multiplication. Then the 


vector {h} is determined from the equation 
[K] {h} = {g} 


by a forward and backward substitution process making use 

Ore the désirablé properties of the LDLT decomposition. Then 
{fi} 1S normalized so as to produce “the next iterant ive .a) 
the reciprocal of the constant of normalization being an 
approximation to we This process is repeated until 
Gonvergence 1S satisfactory. In the computer implementation 
of this procedure, described later, convergence was assumed 
to be satisfactory when two successive such evaluations of 


Z 6 


w, differed by less than uO At this point, the Rayleigh 


quotient was evaluated in the form 
a= (CL)! eva)! Co Cea'eva 7 va! [Mev 


Although the procedure described above could be 
employed with nondiagonal [M], the fact that it was decided 
to diagonalize [M] in order to minimize core storage 
requirements also made the computations involving [M] very 


Simple and economical. 
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IV. PROGRAM ARRANGEMENT 


To facilitate the adaptation of existing mesh generating 
schemes and to incorporate the strong points of existing 
software the general arrangement of the program is modular 
in nature. This also allows future changes to be accomplished 
more easily if desired. 

The actual solution of a problem requires two basic 


steps: mesh generation and frequency evaluation. 


A. MESH GENERATION 

Mesh generation is the process in which the object to 
be analyzed is divided into an assemblage of small elements. 
It is a laborious and tedious task that should not be rushed 
into. Intuition guided by an intelligent study of the 
geometry of the object is often the best approach to uSe in 
choosing a mesh. | 

TRIMEG (Ref. 1, which is discussed in greater detail 
in Appendix C) is a computer program that performs mesh 
generation for three dimensional objects. The program auto- 
matically discretizes a given geometry and upon command will 
punch data decks containing nodal coordinate and element 
connectivity information. The solution routine presented in 
this thesis was written to accept data in the format of 
TRIMEG's output. The use of TRIMEG is not a necessity and any 


Source may be used for the input data. However, the use of 
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TRIMEG to generate data cards for use with NAT1-3D is 
strongly recommended to save effort and to reduce the like- 
lihood of error. The reader is advised to study and 


understand TRIMEG; Ref. 1. 


Bes FREQUENCY EVALUATION 
lies Common Statements 

Throughout the development of this routine a prime 
objective was to make it as simple as possible to use. 
Simplicity can however limit flexibility. To enable the 
user to avoid requesting excessive core space by adjusting 
the core requirements to meet the needs of his particular 
problem, the large data blocks are placed in easily changed 
COMMON storage areas. This adjustment allows the use of 
more nodes for problems with smaller band widths. 

There are four labeled COMMON statements which may 
be used to modify core storage requested by NAT1-3D; this is 
done to avoid requesting more core than necessary and thus 
results in obtaining optimum job priority. The locations 
of the cards and size of the arrays are indicated in the 
listing in Appendix B. The following is a tabulation of the 


Statement names and the subroutines in which they appear. 


ees (Pte nGremOUAD hURM FASTEN, LILD, SOLY 
Sern. INPUMS MERGE 

ER MAIN, NORM 

SOL Weis SO Leal 


2] 








In addition to modifying the COMMON statements the 
dimension of the vector FMODE which appears in the main 
routine should be calculated and introduced in the appropriate 
Meme NSION statement. 

2. Input 

The required data consists of six basic card 
groupings. Their specific formats are listed in Appendix B. 

a. TITLE (one card) 

Any desired problem title may be used. 

b. PROBLEM DATA (one card) 

The card contains the basic problem parameters. 
Here the user chooses the type of shape functions that are 
desired and indicates the number of Gauss points he wishes 
to use in their integration. 

Ce MATERIAL PROPERTIES (one card per material) 

These cards contain the elastic property data 
for the isotropic materials used in the object. The number 
used to represent a material is arbitrary but it must be the 
Same numbering scheme used in the mesh generation. 

d. NODE CARDS (one card per node) 

These cards, which may be punched by TRIMEG, 
contain the node numbers and their global coordinates. 

e. CONNECTIVITY MATRIX (two cards per element) 

These cards, which may be punched by TRIMEG, 
contain the detailed arrangement of all the contributions of 


each element to tne nodal displacements. As the name implies: 


ae 








this matrix connects the elements 


and serves as the control]- 


ling agent which distributes the contribution of each element 


to the global stiffness and mass matrices. 


fie BOUNDARY CONDITION CARDS (one card per 


boundary node) 

The user may impose 
restraining the motion of any node 
direction. The boundary condition 


and the coordinate direction to be 


23 


boundary conditions by 
in any coordinate 
cards indicate the node 


fixed. 





V. EXAMPLE PROBLEMS 


The method usually used to test computer programs is 
to solve a problem that has a Known exact solution to serve 
as a comparison gauge for the computer solution. This 
presents a difficulty as the author is unaware of a truly 
three dimensional classic solution for natural frequency 
determination. Problems A and B, below may seem overly 
simple but they do have known approximate solutions and 
serve to establish a measure of confidence in the integrity 
of the routine. 

The approximate solutions are based on Euler beam theory 
or on the more elaborate Timoshenko beam theory, both of 
which involve replacing the equations of three dimensional 
elasticity (theoretically required to describe the problem) 


by greatly simplified equations. 


BA. CANTILEVER BEAM 


The cantilever beam is a basic design structure which 
may vary in form from the smallest man-made gear tooth to 
the largest trees produced by mother nature. The difficulty 
Oieanalysis also may vary as drastically. The inclusion of 
effects such as shear and rotary inertia significantly 


complicates the frequency equation. 
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The specific structure used in this problem was a beam 
six inches long having a uniform one inch Square cross 
section. Subdivision is indicated in Figure 1. The results 


were as follows: 


Euler Beam Theory w = 5700 rad/sec 
NATI-3D Results w = 5664 rad/sec 
Ratio NAT1-3D/Euler 0.994 


This agrees with Timoshenko theory which predicts a 
frequency correction for shear and rotary inertia of nearly 
unity for this geometry [Ref. 8]. Figure 1 shows the 
convergence of the solution as the number of nodes is 
increased, i.e., the effect of finer subdivision. Note 
however, that the finer discretization involves only the 
"Z" direction. It would be of interest also to investigate 
mamger SUbdiViSi0onNs in the "X° and/or “Y° direction. However, 


this would involve substantial increases in bandwidth. 


B. Spey oUPPORTED BEAM 
The effects of shear deflection and rotary inertia are 
more pronounced in a simply supported beam having the same 


geometry as the previous cantilever. 


Euler Beam Theory Ww 16,000 rad/sec 


NAT1-3D iee2coU rad/sec 


& 
tt 


Timoshenko Beam Theory w = 14,000 rad/sec 
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The Timoshenko solution is based on an assumed shear 
constant value of 0.6. This value is the approximate 


middle of the accepted range of values [Ref. 8]. 


G.. CANTILEVER BEAM WITH COMPLIANT FOUNDATION 

Figure 2 represents a mesh which could be used to 
investigate the effects of root fixity on the natural 
frequency of a cantilever beam having imperfect root fixity. 
The boundary conditions would be applied to the outer wall 
surfaces thus allowing the root to deflect naturally. The 
material properties of the wall and the beam can also be 
varied to measure the effects of dissimilar materials. A 
final solution was not obtainable due to the large core 
requirements arising from even the rather coarse discretiza- 


tion indicated in Figure 2. 
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Discretization of Cantilever to Investigate Root Fixity tffects. 


Figure 2. 





VI. DISCUSS TON 


A. MASS MATRIX 

As noted in chapter II a lumped diagonal mass matrix 
was chosen for use in the solution routine. The choice was 
based upon the reduced core requirements and computational 
advantages inherent to the diagonal matrix. A decision was 
made to form the lumped matrix from the consistent matrix 
because the numerical eva!uation of the consistent mass 
matrix parallels the evaluation of the stiffness matrix 
allowing the two formulations to be coded simultaneously. 
Lumping is accomplished by replacing the diagonal entry in 
each row of the consistent matrix with the sum of the entries 
in that row and zeroing the non-diagonal elements. The 
result could be termed a consistent, lumped matrix. 

Reference 7, previously cited, indicate that this form 
of mass matrix yields excellent results for the lowest 
frequency. It is of some interest at this point to examine 
the point-mass distribution which results from this process 
in the case of an element in the form of a rectangular 
parallelepiped, a case for which physical intuition provides 
Some indication of whether or not the substitution is 


reasonable. 
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Figure 3 Shows the resultant mass distribution for both 
linear and quadratic elements. The distribution for the 
linear element appears logical. The quadratic element 
however is not at all obvious. The allocation of negative 
mass at the element corners would seem to conflict with 
common sense. These results are however, not without 
precedent. Negative corner values result when imposing 
uniform body forces on rectangular elements in a pattern that 
is nearly identical to that shown in Figure 1 of Ref. 6 


foo. 168). 


B. SHAPE FUNCTIONS USED 

In its earlier versions NAT1-3D accommodated quadratic 
shape functions only. This was a compromise between the 
Timitations of Jinear functions in representing complex 
geometry and the large core size requirements of cubic 
elements. A later inclusion of linear shape functions 
resulted in the observation Shee theyawere Coo sSti@t. The 
answers they gave for the fundamental frequency were high 
when the same number of nodes and same nodal arrangement was 
used that corresponded to an accurate solution based on 
quadratic functions. 

This is to be expected since, for the same subdivision 
into elements, the shape functions for quadratic elements 
permit a closer representation of the actual deformation 
pattern than is permitted by using linear elements. In 


effect this corresponds to both types of elements implicitly 
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Figure 3. Allocation of Lumped-Mass Mass Matrix. 


S| 





introducing fictitious and undesirable restraints, which 
increase the calculated frequency above its correct value; 
however, these restraints are greater and have a more serious 
effect in the case of linear elements than they do in the 


case Of quadratic elements. 


ce NUMBER OF GAUSS POINTS 

If N represents the number of Gauss points, the Gaussian 
technique used in this routine will exactly integrate a 
polynomial of degree (2N-1). It has been noted however, in 
calculations employing NAT1-3D that the minimum number of 
Gauss points which can successfully be used with quadratic 
elements is three. When fewer than three Gauss points are 
used it 1s found that accuracy of the frequency evaluation 
is seriously decreased. It is conjectured that this behavior 
is related to the method of introducing the boundary condi- 
tions (by multiplying selected diagonal elements of [BGK] 
by a large factor, "WELD"). The calculations which showed 
the undesirable behavior had WELD = 10'9. It would be 
interesting to look into this matter more fully by reducing 
the value of WELD, by considering alternate methods of 
introducing boundary fixity, and in other ways. Briefly, we 
do not presently have a satisfactory understanding of the 
effect of modifying the number of Gauss points except to 
observe that using NGP = 2 gave unsatisfactory results, while 
using NGP > 3 gave results which were felt to be no more 


pecurate Or reliable than those obtained by using NGP = 3. 
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D. CONVERGENCE 

Our experience has shown that the Stodola iteration 
converges with amazing rapidity, three to four iterations 
usually being sufficient. This can be attributed to the 
fact that the structures considered had lower frequéncies 
which were widely separated from the second frequency. It 
is expected that objects whose lower frequencies are not 
Separated will require more iterations to converge. 

It should be noted that the “initial guess," being a 
unit vector, contains elements of both transverse and axial 
vibrations. If the lowest mode were an axial mode it would 
be reflected in the corresponding elements of the stiffness 
matrix and the iteration should converge to the actual 
lowest frequency be it an axial or transverse mode. Correct 
solutions have been obtained for transverse vibration using 


a pure axial mode as an initial guess. 


an PROGRAM LIMITATIONS 

It could be considered overkill to use this routine on 
a truly two dimensional plane strain or plane stress problem. 
During the development of this routine a 2-D version was 
written but it required exactly the same work on the eerie 
part in the form of data preparation and further required 
the user to recognize that his object could be modeled as a 


plane stress/strain case. Upon considering that an original 
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goal was to limit the assumptionsand idealizations that the 
user has to make, only the three dimensional routine is 
being presented. 

In its present format the program is limited by core 
size. This problem arises when dealing with complex geome- 
tries which require several layers of elements in more than 
one direction. Multidirectional layering results in a 
dramatically rapid growth in the band width of the stiffness 
matrix which is the major core consuming element of the 


entire routine. 
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VII. RECOMMENDATIONS 


A. ELIMINATE STORAGE LIMITATIONS 

The first and most obvious recommendation for further 
improvement of the routine is to remove the core limitations 
imposed on the working area. This can be accomplished by 
the use of direct access devices and a block solving routine 
which sectionalizes the matrices and works in core with only 
one section at a time. 

It should be noted that the implementation of these 
fe@irications will be a difficult task. Substantial coding 
changes will be required and the difficulties encountered 
are anticipated to be orders of magnitude greater than those 


meoctrated with the in-core routine. 


Be PARAMETER STUDY 

It is normally a part of verifying and testing a FEM 
program to perform parametric studies on several problems 
having known exact or approximate results. Time did not 
permit doing such an extensive investigation with NAT1-3D 
and it is recommended that appropriate attention be given to 
a systematic study of the effects of varying the values of 
NGP, WELD, and the degree of the shape functions and of 
varying the fineness of the ultimate discretization in each 
Of the three coordinate directions. The literature on FEM 


describes many such studies and can be used to guide such an 
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Investigation. However, it may be stated that such previous 
investigations have not led to any significant understanding 
of the effects of varying parameters in general. Apparently, 
gemonr yet, it 1S important to carry out such an investigation 


for each new FEM code. 
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APPENDIX A 
PROGRAM DESCRIPTION 


A complete disection and detailed description of the 
Gapgram 1S not practical. This appendix is intended to 
contain information which will assist both a reader who 
desires to pursue a more in-depth study of the problem anda 
a reader who wishes oniy to make an intelligent casual use 


Sietne program. 


A. TERMINOLOGY 
aie tollowing 1S a listing and brief description of 
mportant constant, vector, and matrix names. Functions, 


definitions and dimensions are given if applicable. 


lie Constants 

FACT Value of the Euclidean norm 
of the mode vector, used to 
determine iteration convergence 

NBAND Fol eDanGawaath or the Stadt — 
ness matrix 

NEL Number of elements 

NEQ Number of equations, 
NEQ = 3 x NNODE 

NGP Number of Gauss points used 
in the Gaussian quadrature . 
formula 

NMAT Number of materials 

NNBC Number of bounded nodes 

NNODE Total number of nodes 

NEE Number of nodes per element 
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Ze Vectors (size) 


B Working vectors used to store 

EMODE (NEQ) the mode shape during iteration 

FMODE and while evaluating the 
Rayleigh quotient 


BGM Diagonalized mass matrix stored 
dca VECLOr 

ce Matrices (dimensions) 

BGK(NEQ,NBAND) Gi@bal or system stiffness matix 

COARD(NNODE,3) Matrix containing global 
coordinates of nodes 

ENW(3*NPEL,3*NPEL) Element stiffness or mass matrix 

NCONN(NEL,NPEL+1) Connectivity matrix 

NBC (30,4) Matrix of bounded nodes (i.e., 


nodes on which boundary condi- 
tions are to be imposed) and 
their boundary conditions 


B. PROGRAM ARRANGEMENT 

As indicated in Figure 4 the routine is of modular 
construction. This format allows for maximum ultilization 
of existing software and facilitates future revisions. 

Me MAIN 

The first segments of the MAIN program initiate the 

calculation of the large system matrices. Data is passed to 
the various subroutines via common statements. The iterative 
portion of MAIN is based upon the decomposition of BGK and 
imemrepeated solution of the matrix equation 
meGk] {v} = »° [BGM] {v}. Once convergence has been attained 
MAIN controls the calculation of the Rayleigh quotient which 
serves as both a check and a refinement of the frequency 


obtained by the iteration. 
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ARRANGEMENT OF NAT1-3D 
MAIN 


€-———> | ‘INPUT —> (STOP on empty data file) 
< ? | MERGE | 
hr 
L L v 
QUAD = ELASM aie (2) | 


[Foren | | Form | 
SHAPE FF 
Peer 
| acon | 


<— | FASTEN 


STODOLA ITERATION NAMES in boxes are 
SUD VOULI Ne names 





<— >} LDLT 
C— >» SOLV 
|; ——> | Norm 


RAYLEIGH QUOTIENT 


Figure 4 
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a SUD MDOu tame we UT 

INPUT reads and echo checks the problem data. The 
most common source of errors in a large routine is a simple 
By oodgraphical error in the data deck. INPUT also computes 
the half band width of the system matrices. 

oe Subroutine MERGE [Ref. 4] 

MERGE combines the contributions of each element 
and distributes them to the corresponding global matrix, 
fcr or (BGM. The distribution allocation is controlled 
by the connectivity matrix (NCONN). MERGE serves in a 
sequence controlling capacity throughout the calculation of 
the system stiffness and mass matrices. BGM is computed 
first and during its assembly is stored in the region set 
aside for BGK. Calculation of the weight of the object 
provides a quick check on the mass matrix by giving the user 
an evaluation which he may compare to the actual weight of 
the object. . 

4. Subroutine QUAD [Ref. 4] 

QUAD is divided into two basic sections which 
differ in the form of the shape functions used. Computing 
the mass matrix requires evaluation of the shape functions 
Himee computing the stiffness matrix requires evaluation of 
the derivatives of the shape functions. QUAD selects the 
weighting factors and integration points that correspond to 
the number of Gauss points desired by the user. It then calls 


the appropriate entry of Subroutine FORM which evaluates the 
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integrand of the Gaussian quadrature formula. Next QUAD 
multiplies the result by the Gaussian weighting factor and 
computes the summation. Diagonalization of the mass matrix 
is also accomplished in QUAD. 

Sy SubmolecmeweonncHmlRet. 6) 

ERAStimeronirs the constitutive elasticiwey nmaerix Gor 
each element based on its material which is assumed to be 
ihimear, isotropic and elastic. 

6. Subroutine FORM [Ref. 4] 

FORM evaluates the product of the shape functions 
or the product of their partial derivatives at each selected 
Gauss point in the local coordinate system. Dua! entry 
points are used to distinguish between the computations 
necessary to form the elemental mass matrix and the elemental 
Serrrness matrix. 

ie Subroutine SHAPE [Ref. 4] 

SHAPE performs the actual evaluation of the shape 
functions and their derivatives with respect to the local 
coordinate system. Provisions are included for both linear 
and quadratic shape functions. 

8. Subroutine JACOB [Ref. 4] 
JACOB evaluates the Jacobian matrix as well as its 


inverse and its determinant. 


4] 





oe SUD Rot time FAs) EN 
FASTEN spatially fixes the bounded nodes in the 
desired coordinate direction by multiplying the appropriate 


diagonal elements of [BGK] by oe 


so as effectively to 
Sumoress the associated deflection component. This is a 
commonly used method, in FEM analysis, of introducing fixity 
at boundary nodes. 

Cig Subroutine LDLT [Ref. 5] 

LDLT decomposes the system stiffness matrix [BGK ] 
into the products of three matrices: a lower triangular 
matrix [L] having unit values on the diagonal, a diagonal 
matrix [D], and ry! Ene stranspese 1Ot sev 1 Zz. 

[BGK] = He We eae In the process [BGK] is destroyed and 
replaced by the elements of [D] and the nontrivial elements 
meee. the diagonal matrix is stored in the first column 
and the unit triangular matrix is stored in the remaining 
Space. 

i. Subroutine SOLV [Ref. 5] 

SOLV performs a forward substitution followed by 
amieaek slibstitution to solve the system [BGK] {X} = {B}. 
The answers, {X}, are returned in vector {B}. 

eZ . Subroutine NORM 
NORM calculates the Euclidean norm of the mode 


vector. This norm is equal to the square root of the sum 
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of the squares of the element values in the vector. Each 
element in the vector is then divided by the norm which 


Results in normalizing the vector. 


c CORE ESTIMATION 

The mesh generating scheme chosen by the user should 
inform him of the number of nodes and the bandwidth that 
have resulted from the discretization. TRIMEG clearly lists 
this data. With the aid of the following equation these 
figures may be used to estimate the core requirements which 


result from the mesh. 


a required = 46) = 24 x NNODES 5 -ANBAND eer 


The resultant answer includes the necessary buffers and 
can be used directly as the region specified on the FORTCLG 


control card. 
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AVE ele INI 8 
MESH GENERATION USING TRIMEG 


The purpose of this appendix is to give the user who 
wishes to use TRIMEG a few brief but very important 


Suggestions that could save time and reduce frustration. 


A. BACKGROUND 

TRIMEG [Ref. 1] in an automatic mesh generator for 
use with three dimensional isoparametric elements. It 
1s both quick and efficient and saves the user a lot of 
tedious, time consuming work by punching the bulk of the 


data cards required to solve the vibrational problem. 


B..: COMMENTS 

When using TRIMEG there are a few subtle guides which 
must not be overlooked in order to obtain the most efficient 
aiseretization. 

te: Orientation 

Unless the user has a reasonably good idea of the 

purposes and functions of TRIMEG it is likely that he will 
not use NAT1-3D to the greatest advantage. He is likely to 
orient his super elements or provide for subsequent sub- 
division of super elements into a number of individual 
elements in a diSadvantageous way so as to use up available 
storage while accommodating a smaller total number of elements 


than would have been possible by a more effective arrangement. 
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Accordingly he is urged to study Ref. 2 before undertaking 
Geo GoDlen- we INeparctictiilar, referring to Figtire 13, p. 29, 
Ref. 2, the following general remarks may be of assistance. 

The super elements all have the same orientation given 
by xi, eta, and zeta directions. A super element is divided 
into ROW parts in the eta direction, into COL parts in the 
Mieeadirection, and into SLIGE parts in the zeta direction. 
The body represented by the super elements has, let us say, 
a maximum of nxi super elements in the xi direction, a 
maximum of neta super elements in the eta direction, and a 
maximum of nzeta super elements in the zeta direction. The 
orientation of the super elements should be chosen so that 
the numerical products (ROW)(nxi), (COL)(neta), and (SLICE) 
(nzeta) are in ascending order. This gives minimum bandwidth 
corresponding to the total.number of elements into which the 
body is finally subdivided. 

“Ns Boundary Nodes 

For configurations having curved boundaries it is 

possible and advisable to use 32 potindany nodes to describe 
the super element even when linear or quadratic shape 
functions are to be used. The use of the additional nodes 
allows TRIMEG to more accurately follow the curvature of the 


boundary of the object as it subdivides the super element. 
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NUMBER OF SLICES IN SUPER ELEMENT. 
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